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Chapter  1 


Introduction 

1.1  Motivation 

The  Global  Positioning  System  (GPS)  is  used  for  navigation  and  associated  tasks 
in  a  wide  variety  of  applications.  One  of  the  primary  modes  of  GPS  operation  is  differ¬ 
ential  GPS.  Differential  GPS  allows  more  accurate  positioning  than  stand-alone  GPS. 
Many  new  applications  are  possible  because  of  the  increased  accuracy  obtained 
through  differential  GPS.  For  instance,  a  cruise  missile  operating  with  differential  GPS 
can  be  guided  very  accurately  to  its  target  destination. 

The  benefits  of  differential  GPS  are  apparent,  but  it  must  be  well  understood  in 
order  to  maximize  this  benefit.  There  is  a  need,  therefore,  to  determine  the  positioning 
errors  that  are  associated  with  differential  GPS  operation  for  a  variety  of  cases,  then 
the  user  can  make  knowledgeable  decisions  on  the  implementation  of  the  differential 
technique. 

1.2  GPS  Description 

GPS  signals  are  transmitted  from  the  satellites  at  two  different  frequencies,  LI  at 
1575.42  MHz  and  L2  at  1227.60  MHz.  The  LI  frequency  is  modulated  by  two  pseu¬ 
dorandom  noise  codes,  C/A-code  at  1.023  MHz  and  P-code  at  10.23  MHz.  The  L2  car¬ 
rier  is  only  modulated  by  the  P-code.  Since  the  wavelength  of  the  P-code  is  30  meters 
compared  to  300  meters  for  the  C/A-code,  the  P-code  provides  much  better  accuracy 
for  positioning.  The  P-code  became  known  as  the  Y-code  when  it  was  encrypted  so 
only  authorized  users  could  access  it.  A  P(Y)-code  capable  receiver  must  be  “keyed” 
with  the  proper  information  to  decode  the  incoming  signals  and  obtain  P-code  accu¬ 
racy.  [7] 
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For  C/A-code  users,  the  performance  is  degraded  due  to  a  purposeful  degradation 
in  the  quality  of  broadcast  orbits  and  satellite  clock  dithering.  This  degradation  is 
called  selective  availability  (SA).  [24]  An  advantage  of  P-code  operation  is  that  per¬ 
formance  is  not  affected  by  selective  availability. 

Positioning  is  accomplished  by  determining  the  time  it  takes  for  the  GPS  signals 
from  each  satellite  to  reach  the  receiver.  This  time  provides  an  estimate  of  the  distance 
from  satellite  to  receiver,  and  “triangulating”  using  the  signals  from  four  or  more  satel¬ 
lites  provides  the  receiver’s  position  in  three  dimensions.  The  transmit  and  receive 
times  are  corrupted  by  clock  errors,  so  a  true  range  measurement  is  not  possible. 
Therefore,  the  distance  measurement  from  the  receiver  to  each  satellite  is  called  the 
pseudorange. 

The  “triangulation”  is  not  perfect,  because  the  satellites  are  not  located  optimally 
with  respect  to  the  receiver.  Additional  error  is  contained  in  the  navigation  solution 
due  to  the  position  of  the  satellites.  The  instantaneous  measure  of  this  error  is  called 
the  Geometric  Dilution  of  Precision  (GDOP).  It  represents  the  rms  error  realized  in  the 
navigation  solution  for  unit  error  in  the  triangulation  measurements.  GDOP  accounts 
for  dilution  of  precision  in  three  dimensions  and  time.  GDOP  can  be  broken  down  into 
the  following  components:  vertical  (VDOP),  horizontal  (HDOP),  time  (TDOP),  posi¬ 
tion  (PDOP).  PDOP  reflects  the  dilution  of  precision  in  three  dimensions.  [20]  This 
can  be  further  decomposed  into  XDOP,  YDOP,  and  ZDOR 

Incoming  GPS  signals  are  delayed  when  they  pass  through  the  earth’s  atmosphere. 
The  ionosphere  is  one  of  the  primary  causes  for  this  delay.  A  P-code  receiver  can  track 
both  the  LI  and  L2  signals.  Since  these  signals  are  at  different  frequencies,  the  iono¬ 
spheric  delays  will  be  different  for  each.  An  estimate  for  the  ionospheric  delay  can  be 
obtained  utilizing  the  difference  between  the  LI  and  L2  delays.  Thus  a  P-code  receiver 


14 


provides  a  dual-frequency  ionospheric  correction  term  that  is  subtracted  from  pseudo¬ 
ranges.  Since  the  LI  and  L2  delays  are  obtained  from  tracking  loops  on  the  signals, 
additional  noise  enters  the  navigation  solution  when  the  dual-frequency  correction  is 
used.  [10] 

GPS  satellites  broadcast  almanac  files  that  contain  the  orbit  parameters  required  to 
predict  the  satellites’  positions.  These  files  are  updated  every  week  in  order  to  main¬ 
tain  the  accuracy.  All  GPS  data  is  referred  to  by  week.  The  system  clock  is  reset  to 
zero  every  week  on  midnight  Saturday,  so  the  correct  date  is  found  by  referencing  the 
GPS  week  in  the  data  files. 

1.3  Differential  GPS 

As  shown  in  Figure  1.1,  positioning  using  differential  GPS  (DGPS)  involves  the 
use  of  two  GPS  receivers.  One  receiver  is  located  at  a  fixed  position  that  is  accurately 
surveyed,  while  the  other  is  allowed  to  move  around.  The  goal  is  to  accurately  deter¬ 
mine  the  position  of  the  “roaming”  receiver. 
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Since  the  position  of  receiver  A  is  surveyed,  corrections  can  be  determined  for  the 
pseudoranges  that  generate  the  GPS  position.  These  corrections  are  then  sent  to 
receiver  B.  The  roaming  receiver  uses  the  corrections  from  receiver  A  to  update  or 
correct  its  position.  This  procedure  is  highly  accurate  with  small  baseline  lengths. 
Errors  increase  as  the  baseline  increases  or  as  a  time  delay  is  introduced  between  the 
measurement  and  use  of  the  corrections.  Such  a  time  delay  can  be  the  result  of  the 
desire  for  covert  operation. 

Using  DGPS  can  help  to  eliminate  errors  due  to  the  satellite  clock  drift,  ephemeris 
errors,  and  propagation  delays.  Errors  remaining  include  receiver  interchannel  biases, 
receiver  noise  and  quantization  errors,  multipath,  residual  propagation  delays,  and 
errors  due  to  time  delays  in  the  pseudorange  corrections.  [4] 

1.4  Previous  Work 

Several  sources  were  found  that  discussed  DGPS  errors.  Most  of  the  previous 
work  was  performed  with  C/A-code  receivers.  Error  budgets  are  cited  in  [4,  15,  17] 
and  depend  on  the  type  of  receiver  being  used.  For  the  most  part,  these  budgets  were 
the  specifications  for  the  devices  or  analytic  estimates  which  do  not  always  correlate 
with  actual  performance.  The  next  section  provides  one  of  these  error  budgets.  An 
analytic  development  of  the  individual  error  sources  was  provided  in  [9]. 

Two  references  were  found  that  model  GPS  and  DGPS  errors  as  first  order  Markov 
processes.  [7,  23]  A  Markov  model  was  chosen  because  many  physical  phenomena 
exhibit  behavior  that  fit  such  models  reasonably  well.  No  GPS  data  was  used  in  deriv¬ 
ing  these  models. 

One  reference  used  simulated  data  because  receiver  data  for  the  desired  scenarios 
was  unavailable.  [14]  Some  studies  used  C/A-code  data  for  their  analysis,  including 
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data  found  on  the  Internet.  [5,  13]  One  study  did  use  P-code  receivers,  but  this  study 
examined  the  DGPS  errors  for  a  specific  geometry  and  did  not  develop  a  model.  [25] 
This  thesis  utilizes  data  segments  with  common  satellites  and  processes  them  inde¬ 
pendently.  This  prevents  jumps  in  the  solutions  due  to  satellite  switches  from  affecting 
the  results.  A  statistical  model  for  the  DGPS  errors  is  provided. 

1.5  Error  Budgets  and  Modeling 

Errors  for  a  P-code  receiver  are  typically  cited  at  about  10  meters  RMS  in  each 

axis.  Table  1.1  shows  the  error  budgets  for  static  GPS  operation,  and  Table  1.2  shows 
an  error  budget  for  Differential  GPS.  These  budgets  are  specifications  for  performance 
for  co-altitude  receivers  at  a  separation  of  50  nmi  (92.6  km).  They  were  published  by 
the  Air  Force  GPS  Range  Applications  Joint  Program  Office.  [4,  15] 

Table  1.1:  GPS  Static  Error  Budget 


Error  Source 

C/A-Code  Predicted  Error 
(meters) 

P-Code  Predicted  Error 
(meters) 

Satellite  Clock  Error 

3.05 

3.05 

Ephemeris  Error 

2.62 

2.62 

Ionospheric  Delay  Error 

6.40 

0.40 

Tropospheric  Delay  Error 

.40 

0.40 

Receiver  Noise 

2.44 

0.24 

Receiver  Interchannel  Bias 

0.61 

0.15 

Multipath 

3.05 

1.22 

UERE  (RMS) 

8.54 

4.25 

RMS  Horizontal  Position 
Error  (HDOP=  1.5) 

12.81 

6.37 

RMS  Vertical  Position 
Error  (VDOP  =  2.5) 

21.34 

10.62 

Total  RMS  Error 

24.89 

12.38 
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Table  1.2:  Differential  GPS  Error  Budget 


Error  Source 

C/A-Code  Predicted  Error 
(meters) 

P-Code  Predicted  Error 
(meters) 

Satellite  Clock  Error 

0 

0 

Ephemeris  Error 

0 

0 

Residual  lonospheric/Tro- 
pospheric  Delay  Error 

0.15 

0.15 

Receiver  Noise 

2.44 

0.24 

Receiver  Interchannel  Bias 

0.61 

0.15 

Multipath 

3.05 

1.22 

UERE  (RMS) 

3.95 

1.26 

RMS  Horizontal  Position 
Error  (HDOP  =  1.5) 

5.93 

1.89 

RMS  Vertical  Position 
Error  (VDOP  =  2.5) 

9.88 

3.15 

Total  RMS  Error 

11.53 

3.68 

The  DGPS  errors  will  change  as  a  function  of  baseline  distance  and  the  time  delay 
in  the  correction  inputs.  Previous  work  has  estimated  the  change  with  baseline  dis¬ 
tance  to  be  a  linear  relationship.  [13]  One  of  the  major  challenges  in  this  project  will 
be  estimating  long  term  behavior  from  many  segments  of  short  term  data. 


1.6  GPS  Data 

Several  data  sets  were  utilized  in  this  thesis.  Draper  Laboratory  collected  all  the 
data.  Three  data  collections  were  performed:  East  Coast  experiments,  West  Coast 
experiments,  and  zero  baseline,  rooftop  experiments.  The  first  two  experiments  were 
performed  under  a  previous  effort.  They  provide  GPS  data  at  locations  with  varying 
separations  so  several  baselines  are  available.  The  experiment  names  refer  to  the 
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region  where  the  experiments  were  performed,  and  the  objectives  for  both  were  the 
same.  The  rooftop  experiment  was  performed  as  part  of  this  thesis.  The  primary  bene¬ 
fit  was  to  be  a  characterization  of  the  receiver  noise.  An  unexpected  benefit  was  char¬ 
acterization  of  errors  due  to  multipath.  This  is  possible  because  the  baseline  is  zero, 
and  both  receivers  take  inputs  from  the  same  antenna.  These  experiments  will  be 
described  in  more  detail  in  Chapter  2.  The  results  from  all  three  experiments  will  be 
combined  to  form  a  DGPS  model  from  many  different  cases. 

1.7  Thesis  Objectives 

This  thesis  will  examine  the  errors  in  differential  GPS  as  a  function  of  baseline  dis¬ 
tance  between  the  two  receivers  and  the  time  delay  between  when  the  corrections  are 
computed  and  when  they  are  used.  The  final  results  will  include  data  from  P-code 
receivers  only.  A  model  for  the  errors  will  be  developed  and  analyzed  to  determine 
how  well  it  fits  the  errors  observed  in  actual  differential  GPS  experiments. 

The  thesis  is  organized  as  follows.  Chapter  2  discusses  the  experiments  and  the 
data  obtained.  Chapter  3  presents  methods  used  in  the  processing  of  the  data  and  mod¬ 
eling  of  the  errors.  In  Chapter  4,  the  results  from  each  experiment  is  shown  including 
the  performance  of  the  models  developed.  Chapter  5  is  a  summary  of  the  thesis  and 
conclusions  drawn  from  the  results.  Suggestions  for  future  study  are  also  discussed. 
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Chapter  2 


Data  Collection  and  Formats 

This  chapter  describes  the  data  that  was  used  to  develop  the  model  for  DGPS 
errors.  All  of  this  data  was  collected  by  Draper  Laboratory.  Much  of  it  comes  from 
keyed  P-code  receivers  that  is  unaffected  by  selective  availability. 

2.1  Draper  Rooftop  Experiments 

2.1.1  Experiment  Description 

The  goal  of  the  rooftop  experiments  was  to  obtain  data  with  a  zero  baseline  and  to 
get  an  estimate  of  receiver  noise.  A  GPS  antenna  was  placed  on  a  tripod  on  the  rooftop 
of  Draper.  Two  Trimble  TANS  receivers  in  a  laboratory  on  the  4th  floor  received  the 
inputs  from  the  single  antenna.  Since  the  receivers  were  getting  the  exact  same  inputs, 
the  only  difference  in  their  solutions  was  due  to  the  receivers’  internal  processing.  The 
baseline  is  zero  meters,  and  the  time  delay  errors  were  examined  independently  of 
basehne  dependent  errors.  Data  was  recorded  every  10  seconds. 

The  TANS  receivers  are  keyed  P-code  receivers  which  are  resistant  to  selective 
availability.  The  position  errors  are  expected  to  be  well  within  the  specified  value  of  16 
meters.  A  bug  in  the  translation  program  written  for  the  TANS  receivers  caused  errors 
in  the  recorded  ephemeris  data.  A  Novatel  receiver  was  positioned  on  the  rooftop  to 
collect  data  for  the  same  periods  as  the  TANS  receivers.  The  Novatel  translation  pro¬ 
gram  worked,  so  it  was  used  as  the  source  of  broadcast  ephemeris  data. 
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2.1.2  Data  Description 

The  first  set  of  data  was  taken  on  Monday,  November  14th,  1995,  23:44:50  to 
Tuesday,  November  15th,  22:33:10.  This  was  GPS  week  775.  This  set  contains 
approximately  24  hours  of  data.  Unfortunately,  one  of  the  computers  controlling  the 
TANS  receivers  failed  to  store  the  collected  data.  All  data  from  one  receiver  was  lost. 
It  was  thus  not  possible  to  estimate  receiver  noise  from  this  data.  However,  the  data 
from  the  remaining  receiver  is  very  useful  for  the  time  delay  error  analysis,  because 
the  data  can  be  copied  and  shifted  in  time.  In  many  trials  comparing  receiver  A  vs. 
receiver  B  (two  different  receivers)  with  receiver  A  vs.  receiver  A  (the  same  receiver), 
the  time  shifted  errors  are  statistically  identical.  Using  the  same  receiver  for  the  time 
delay  study  eliminates  any  sample  mean  in  the  data. 

The  second  set  of  data  was  taken  on  Thursday,  November  17th,  1995,  17:02:38  to 
22:13:40.  This  was  also  GPS  week  775.  This  set  contains  approximately  5  hours  of 
data.  This  time  both  computers  stored  the  data,  so  comparisons  between  the  two 
receivers  could  be  performed.  Unfortunately,  a  gap  of  about  1.5  hours  occurred  in  the 
middle  of  the  data.  The  gap  occurred  when  neither  receiver  was  able  to  maintain  signal 
lock  on  more  than  three  satellites.  The  reason  for  the  gap  has  not  been  determined. 
Possible  explanations  include  objects  blocking  the  line  of  sight  to  the  satellite,  an  inac¬ 
tive  GPS  satellite,  or  a  satellite  that  was  not  correctly  transmitting  at  that  time.  It  is 
unlikely,  however,  that  any  nearby  structure  could  have  blocked  the  satellites  for  1.5 
hours. 

The  receivers  were  configured  to  track  the  same  four  satellites.  This  was  done  so 
the  position  solution  from  both  receivers  would  have  common  error  sources.  Using 
solutions  from  different  satellites  would  mask  the  differences  we  are  trying  to  observe. 
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The  receivers  were  also  configured  to  hold  a  set  of  four  satellites  until  the  GDOP 
became  greater  than  6.  At  that  time,  both  receivers  switched  to  a  new  set  of  four  satel¬ 
lites. 

2.2  East  Coast  Experiments 
2.2.1  Experiment  Description 

The  East  Coast  experiments  were  conducted  in  December,  1993.  These  experi¬ 
ments  included  both  C/A  and  P-code  receivers.  Table  2.1  shows  the  location  and  type 
of  the  receivers  used.  Data  collected  from  Ft.  Stewart,  Georgia  was  not  used  due  to 

Table  2.1  East  Coast  Receiver  Locations  and  Types 


Receiver  Location 

Receiver 

Type 

Cambridge,  MA 

P-code 

Portsmouth,  NH 

P-code 

Langley  AFB,  VA 

P-code 

Willow  Run,  MI 

C/A-code 

problems  with  the  data  collection.  The  East  Coast  baselines  are  shown  in  Table  2.2. 

Table  2.2  East  Coast  Baselines 


Receiver  A 
Location 

Receiver  B 
Location 

Baseline  Length 
(km) 

Cambridge,  MA 

Portsmouth,  NH 

84 

Langley  AFB,  VA 

Cambridge,  MA 

740 

Langley  AFB,  VA 

Portsmouth,  NH 

820 

Langley  AFB,  VA 

Willow  Run,  MI 

838 

Willow  Run,  MI 

Cambridge,  MA 

1022 

Willow  Run,  MI 

Portsmouth,  NH 

1052 
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The  C/A-code  receiver  at  Willow  Run  was  a  Novatel  95 IR.  The  selective  availability 
in  the  data  from  this  receiver  was  removed  by  Aerospace  Corporation.  The  P-code 
receivers  were  Trimble  4000SSE  unencrypted  receivers.  The  P-code  receivers  were 
able  to  track  P-code  even  though  unencrypted,  because  P(Y)-code  encryption  was  not 
yet  in  effect. 

2.2.2  Data  Description 

Seventeen  different  sessions  of  data  were  collected.  All  the  receivers  were  oper¬ 
ated  at  the  same  time.  Some  of  the  sessions  were  lost  due  to  data  collection  problems. 
The  data  was  searched  for  segments  where  all  receivers  had  at  least  twenty-five  min¬ 
utes  of  time  tracking  the  same  four  satellites.  Five  such  segments  were  found  in  ses¬ 
sions  8,  10,  14,  and  16.  These  five  data  segments  were  used  to  form  time  delay 
differences  as  well  as  to  estimate  the  differences  due  to  the  six  baseline  separations. 

2.3  West  Coast  Experiments 
2.3.1  Experiment  Description 

The  West  Coast  experiments  were  conducted  in  March,  1994.  Keyed  P-code  Trim¬ 
ble  4000SSE  receivers  were  used  at  all  locations.  Only  two  such  receivers  were  avail¬ 
able,  however,  so  each  baseline  had  to  be  instrumented  separately.  Data  was  collected 
from  four  baselines  as  shown  in  Table  2.3. 


Table  2.3  West  Coast  Baselines 


Receiver  A 
Location 

Receiver  B 
Location 

Baseline  Length 
(km) 

San  Diego,  CA 

Long  Beach,  CA 

149 

Stockton,  CA 

Fallon,  NV 

282 

Ft.  Irwin,  CA 

Fallon,  NV 

496 

San  Diego,  CA 

Fallon,  NV 

758 
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2.3.2  Data  Description 

Sixteen  sessions  of  data  were  collected  on  the  west  coast,  including  four  sessions 
per  baseline.  Data  collection  problems  resulted  in  the  loss  of  two  sessions.  For  the 
remaining  sessions,  data  segments  of  25  minutes  or  greater  with  the  two  receivers 
tracking  the  same  four  satellites  were  selected.  Each  session  had  from  1  to  5  such  seg¬ 
ments.  These  segments  were  again  used  to  form  the  time  delay  differences  and  to  esti¬ 
mate  the  error  due  to  the  four  baseline  separations. 

2.4  Data  Anomalies 

Several  interesting  qualities  or  anomalies  have  been  observed  in  the  collected  data 
sets.  These  are  documented  to  stimulate  future  discussion  and  investigation.  The  gap 
in  the  5-hour,  zero  baseline  rooftop  data  has  already  been  mentioned.  From  discus¬ 
sions  with  experienced  GPS  users,  satellite  outages  have  occurred  on  occasion. 

The  GPS  navigation  solution  jumps  at  every  satellite  switch.  This  is  an  expected 
attribute,  since  it  is  a  function  of  all  four  satellites  being  tracked.  One  of  the  primary 
goals  of  the  data  reduction  was  to  deal  appropriately  with  these  jumps.  Figure  2.1 
shows  the  data  from  a  zero  baseline  5-hour  rooftop  session.  The  graphs  show  the  com¬ 
ponents  of  the  difference  between  one  receiver’s  solution  and  presumed  truth.  The  ver¬ 
tical  lines  show  the  places  where  the  satellite  switches  occur.  Notice  the  jumps  in  the 
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Figure  2.2  shows  a  close-up  of  the  last  satellite  break  in  the  zero  baseline  data. 
Note  the  fairly  large  jumps  in  both  x  and  z  directions. 


GPS  Time  (hours) 

Figure  2.2  Satellite  Break,  Zero  Baseline  Data 

Figure  2.3  shows  a  close-up  of  another  session  of  zero  baseline  data.  The  vertical 
lines  here  are  not  the  satellite  switches.  They  mark  each  measurement  made  by  the 
GPS  receiver.  These  occur  every  10  seconds.  A  satellite  switch  occurred  at  1 14.32  sec¬ 
onds.  Note  the  data  gap  before  the  switch.  This  occurred  several  times  in  the  collected 
data,  and  is  presumably  due  to  a  delay  in  locking  onto  the  new  satellite. 

The  errors  often  exhibited  cyclic  behavior.  A  change  in  the  frequency  and  magni¬ 
tude  of  these  oscillations  is  observable  at  the  satellite  switches.  The  oscillations  will  be 
primarily  attributed  to  multipath.  The  changes  in  the  oscillations  make  sense  since  the 
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new  satellite  set  provides  a  new  multipath  geometry.  This  will  be  discussed  further  in 
the  next  chapter. 
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Figure  2.3  Satellite  Switch  and  Data  Gap,  Zero  Baseline  Data 


Figure  2.4  shows  a  sample  of  the  West  Coast  data.  Again  the  vertical  lines  mark 
each  satellite  switch.  Jumps  in  one  or  more  components  are  apparent  at  each  satellite 
switch. 


Figure  2.5  shows  a  close-up  of  the  West  Coast  data.  The  vertical  lines  mark  each 
GPS  measurement.  Data  was  taken  every  30  seconds.  A  satellite  switch  occurred  at  the 
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jump  in  the  middle  of  the  data  set  at  63.25  seconds.  Note  once  again  the  jump  in  the 
solution  that  occurred  at  the  satellite  switch. 
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Figure  2.5  Close-up  of  West  Coast  Data 


Figure  2.6  shows  an  example  with  only  the  y  direction  plotted.  Note  that  satellite 
switches  occurred  at  86.8  and  87.1  seconds.  The  changes  in  amplitude  and  frequency 
of  the  oscillations  in  the  navigation  solution  are  more  observable  in  this  case.  Analysis 
will  show  that  primary  oscillation  periods  change  from  71.4  seconds  to  90.9  seconds 
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to  125  seconds  for  the  satellite  switches  in  this  example.  This  was  determined  by 
applying  a  fast  Fourier  transform  to  the  data. 


The  GPS  navigation  solution  was  observed  to  jump  at  points  other  than  satellite 
switches.  In  one  case,  a  jump  in  z  position  of  7  meters,  y  of  3  meters,  and  x  of  0.75 
meters  was  observed.  Such  jumps  were  rare,  and  the  cause  is  unknown.  One  explana¬ 
tion  is  that  the  satellite  ephemeris  was  updated  at  those  points.  Figure  2.7  shows  a 
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close-up  of  a  jump  in  the  solution  when  there  was  no  satellite  switch.  The  reader  may 
have  noticed  this  jump  in  Figure  2.1,  towards  the  very  end  of  the  data  set. 


GPS  Time  (hours) 

Figure  2.7  Navigation  Solution  Jump  Without  Satellite  Switch 

Another  anomaly  of  note  was  an  occasional  large  and  sometimes  persistent  differ¬ 
ence  between  the  position  solution  from  the  two  receivers.  Given  that  the  receivers 
were  being  driven  by  the  same  antenna,  the  difference  had  to  have  arisen  from  some 
internal  difference  between  the  two  receivers.  Sometimes  these  differences  only  lasted 
for  one  measurement,  sometimes  they  lasted  for  100  seconds  (10  measurements).  A 
GPS  receiver  is  a  complex  non-linear  device.  Such  behavior  is  not  surprising  and 
could  arise  for  a  number  of  reasons.  Included  are  low  S/N  for  a  particular  satellite  or 
some  receiver  moding  change.  The  level  of  process  noise  being  introduced  into  the 
receiver  Kalman  filter  was  such  that  they  were  essentially  delivering  point  solutions. 
This  was  done,  of  course,  so  that  no  “artificial”  time  correlation  would  be  introduced 
into  the  data.  The  persistence  of  the  anomalies  cannot  thus  be  blamed  on  the  naviga- 
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tion  filter.  It  should  be  remembered  that  the  receivers  were  being  forced  to  hold  onto  a 
set  of  satellites  until  the  GDOP  had  risen  to  7.  This  would  tend  to  make  the  solution 
more  sensitive  to  pseudo-range  errors  than  would  normally  be  the  case. 

The  ephemeris  data  collected  by  one  of  the  Trimble  receivers  was  missing  some  of 
the  satellite  information.  It  appeared  to  start  taking  data  correctly,  but  then  ephemeris 
records  would  show  only  information  from  about  three  of  the  satellites.  Fortunately, 
the  ephemeris  from  the  Novatel  receiver  was  available,  so  this  loss  did  not  affect  the 
results.  However,  no  reason  for  the  data  loss  has  been  determined.  All  other  elements 
of  the  data  from  that  receiver  seemed  fine. 

The  Novatel  receiver  ephemeris  changed  GPS  week  from  775  to  776.  This 
occurred  on  Thursday,  November  17th,  at  6:00  p.m.  The  change  in  GPS  week  was 
expected  at  midnight  on  Saturday.  The  results  looked  fine  and  confirmed  that  the 
change  had  actually  been  made.  The  almanac  data  from  week  776  was  used,  and  the 
times  from  -55  to  -49  hours  corresponded  to  113  to  119  hours  of  GPS  week  775.  It  is 
noted  that  in  the  Novatel  ephemeris  file,  the  almanac  entries  read  week  776,  while  the 
navigation  solution  entries  still  read  week  775. 
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Chapter  3 


Data  Processing 

3.1  Introduction 

The  goal  is  to  observe  and  model  the  time  and  distance  correlated  behavior  of 
DGPS  errors.  One  such  model  is  a  first  order  Markov  process, 

( 

MSE  =  a^ps  1-e^  (3.1) 

V  y 

where  MSE  is  the  mean  square  error,  o^pg  is  the  variance  of  the  Markov  process,  x  is 

the  time  correlation  constant,  and  x  is  the  distance  correlation  constant.  A  notional 

c 

plot  of  this  model  is  shown  in  Figure  3.1.  We  will  focus  for  the  most  part  on  Markov 
models  because  of  their  utility  in  navigation  error  analysis.  The  basic  underlying  phe¬ 
nomena  which  cause  these  errors  may  introduce  correlation  between  the  time  and  dis¬ 
tance  behavior.  An  example  would  be  a  traveling  ionospheric  disturbance.  Over  an 
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ensemble  of  cases,  however,  the  correlation  is  expected  to  be  zero.  It  will  be  so 
assumed  for  this  study. 
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Figure  3.1  Notional  Mean  Square  Error  Model  for  DGPS 

3.2  Inputs  to  the  Processing 

The  first  step  in  the  processing  was  to  remove  corrections  for  ionospheric  delay. 
This  eliminates  the  noise  from  the  L2  channel,  reducing  the  noise  level  by  a  factor  of 
^ .  Differencing  outputs  from  the  two  receivers  will  cancel  the  ionospheric  errors  to 
the  extent  that  they  are  the  same.  This  is  a  common  practice  in  DGPS  operation.  Fur¬ 
ther  characterization  of  this  residual  ionospheric  error  is  one  of  the  goals  of  the  study. 
Removal  of  the  ionospheric  correction  was  actually  done  prior  to  the  start  of  this  study. 
[10]  The  format  for  the  data  is  shown  in  Table  3.1. 

Table  3.1  GPS  Data  File  Format 


1 

2 

3 

4 

5 

6 

7 

8 

9 

Time 

X 

Y 

Z 

#  Sats 

Pml 

Prn2 

Pm3 

Prn4 

10 

11 

12 

13 

14 

15 

16 

17 

18 

Pm5 

Pm6 

GDOP 

PDOP 

HDOP 

VDOP 

TDOP 

Status 

Clock 

The  measurements  from  the  two  receivers  A  and  B  were  edited  by  deleting  mea¬ 
surements  for  which  the  time  tags  could  not  be  matched.  In  addition,  measurements 
were  used  only  when  the  receivers  indicated  that  the  solution  was  good  (status  flag  = 
1). 

Satellite  switch  times  were  computed  in  order  to  split  the  data  into  different  sets 
that  used  the  same  four  satellites  for  the  entire  set.  This  intentionally  eliminates  satel¬ 
lite  switching  as  one  on  the  sources  of  decorrelation.  A  negative  aspect  of  this  strategy 
is  that  we  will  be  trying  to  estimate  long  term  behavior  from  a  number  of  short  term 
data  segments.  This  will  be  discussed  in  more  detail  later.  The  data  segments  were 
samples  from  what  is  presumed  to  be  a  zero  mean  process.  The  small  sample  size 
introduced  a  sample  mean  into  the  data.  This  mean  was  removed  before  determining 
the  correlation  between  different  data  sets. 

Another  strategy  for  getting  as  close  to  the  underlying  phenomenon  as  possible, 
i.e.  errors  in  the  pseudo-ranges,  was  to  normalize  the  errors  in  the  position  solution 
components  using  the  dilution  of  precision  in  each  component,  XDOP,  YDOP,  and 
ZDOR  These  values  were  computed  by  utilizing  a  GPS  constellation  simulation  and 
inputting  the  measurement  times  for  each  set  of  data  and  the  satellites  tracked  by  each 
receiver.  [6]  The  program  makes  use  of  the  almanac  file  taken  from  the  Novatel 
receivers.  The  DOP  values  were  obtained  from  the  diagonal  of  the  (h’^H)  matrix. 
The  resulting  DOPs  were  used  to  normalize  the  x,  y,  and  z  errors. 

For  the  rooftop,  zero  baseline  data  the  GPS  antenna  was  not  placed  at  the  surveyed 
marker.  (The  marker  is  in  the  “shadow”  of  a  large  vertical  metal  wall.)  The  placement 
was  in  an  attempt  to  minimize  multipath  effects  on  the  solution.  Therefore,  the  truth 
value  for  the  position  (x,  y,  z)  was  defined  to  be  the  average  position  for  all  the  data 
sets.  This  differed  from  the  surveyed  position  by  5.6590  meters. 
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Ysurv  +  5.8299  meters,  and  +  18.7744  meters.  These  values  were  in  the  right 
direction  and  distance  from  the  marker  -  as  well  as  could  be  determined  by  eyesight. 
No  effort  was  made  to  actually  measure  the  distance,  since  it  was  noted  that  the  survey 
marker  position  was  itself  only  accurate  to  about  a  meter.  [21] 

3.3  Definition  of  the  Error  and  Difference  in  Errors 

This  section  defines  the  errors  that  were  utilized  in  the  autocorrelation  processing. 

The  quantity  of  interest  is,  finally,  the  difference  in  error  between  the  time  shifted 
errors  from  receiver  A  and/or  B. 

^diff  =  (^err  +  At)  -  (t)  )  /XDOP  (3.2) 

The  RSS  difference  is 

“  V^diff  ydiff  ^diff 

The  individual  components  of  the  error  at  each  receiver  are  defined  by: 

Xerr  =  (3-4) 

These  individual  errors  are  split  into  six  sets  where  all  the  points  in  a  set  use  the 
exact  same  four  satellites.  The  RSS  for  each  set  was  found  using  Equation  (3.3). 

3.4  Auto  and  Cross  Correlation 

The  first  strategy  for  characterizing  the  behavior  of  the  differences  was  to  examine 
the  correlation  in  the  rooftop  data.  The  cross  correlation  and  the  autocorrelation  of  the 
rooftop,  zero  baseline,  differences  was  computed.  Since  the  two  receivers  were  being 
driven  from  the  same  antenna,  and  the  receivers  were  only  contributing  white  noise, 
the  cross  and  auto  correlation  were  expected  to  both  be  measures  of  the  same  thing, 
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the  behavior  of  the  errors  in  pseudorange  with  time.  This  analysis  was  done  for  the 
component  differences  and  the  RSS  difference  for  each  of  6  data  sets. 

The  cross  correlation  between  receiver  A  and  B  differences  was  found.  The  sample 
mean  was  subtracted  from  the  data  before  the  correlation  algorithm  was  applied.  Cor¬ 
relations  were  found  for  x,  y,  z,  and  the  RSS  differences  for  each  set  from  1  to  6  using 
Equation  (3.2)  and  (3.3). 

The  sigma  and  tau  values  were  estimated  for  each  set.  The  variance,  a  ,  is  the  cen¬ 
ter  value,  since  the  inputs  were  of  zero  mean.  Tau  is  the  value  on  the  x  axis  that  corre- 

2 

sponds  to  a  y  value  of  — . 

e 

A  Markov  process  has  a  probability  distribution  that  is  dependent  only  on  the 
value  at  one  point  immediately  in  the  past.  The  process  can  be  described  by 
Equation  (3.5),  where  w  is  white  noise  and  p  is  the  reciprocal  time  constant.  [12] 

/Vv 

^  +  p  (t)  X  =  w  (3.5) 

The  autocorrelation  of  a  zero  mean  Markov  process  is  given  by  Equation  (3.6). 
The  correlation  time,  or  1/e  point,  is  found  at  1/p .  [12] 


(.)  = 


A  single  variance  and  time  constant  for  the  six  set  ensemble  was  found  in  two 
ways:  a)  by  finding  the  correlation  of  the  weighted  average  of  the  data  in  the  six  sets, 
and  b)  taking  the  weighted  average  of  the  variance  and  time  constant  found  from  the 
correlation  of  each  of  the  six  sets. 


The  results  of  the  cross-correlation  process  were  disappointing  for  two  reasons. 
The  first  is  that  a  tremendous  amount  of  data  is  required  to  get  accurate  results  using 
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this  technique.  Brown  and  Hwang  [7]  present  the  limit  on  experimental  autocorrela¬ 
tion  as 


(3.7) 


where  i  is  the  time  constant,  V  (x)  is  the  experimentally  determined  autocorrelation 

P  X 

2 

function,  T  is  the  time  length  of  the  experimental  record,  and  a  is  the  variance  of  the 
process. 

For  a  requirement  of  10%  accuracy  on  the  estimation  of  the  autocorrelation  func¬ 
tion,  the  experimental  autocorrelation  should  have  a  standard  deviation  less  than  0. 1 


times  a  ,  or  0.01  times  a  .  Therefore, 


Var 


=  0.01 ,  and  the  limit  becomes 


T  = 


0.01  p 


(3.8) 


for  an  accuracy  of  10%.  If  the  time  constant  is  1800  seconds,  T  =  360,000  sec¬ 
onds.  In  other  words,  4.2  days  of  data  is  required.  In  general,  the  data  required  for  10% 
accuracy  equals  200  times  the  time  constant.  A  time  constant  of  250  seconds  requires 
almost  14  hours  of  data  for  10%  accuracy.  The  development  of  the  accuracy  limit  is 
presented  in  Appendix  B. 

Figure  3.2  is  a  graphic  expression  of  Equation  (3.7).  It  shows  the  plot  of  amount  of 
data  vs.  percent  error  in  the  autocorrelation  function  for  a  process  with  a  time  constant 
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of  1800  seconds.  Again,  note  that  about  four  days  worth  of  data  is  required  in  order  to 


Amount  of  Data  vs  Percent  Error  for  Time  Constant  of  1800  Seconds 


obtain  about  10%  error  on  the  autocorrelation  function.  For  our  problem,  we  must 
have  data  where  the  satellites  are  all  the  same  for  each  set  of  data.  This  limits  the  data 
sets  to  about  one  hour  before  a  satellite  switch  occurs.  The  second  reason  for  abandon¬ 
ing  correlation  as  the  analysis  method  was  that  it  does  not  easily  allow  the  resolution 
of  multiple  Markov  processes.  It  will  be  found  that  there  were  indeed  two  separable 
Markov  processes  with,  obviously,  different  time  constants. 

3.5  Time  Delay  Error  Analysis 

An  alternative  method  for  observing  and  characterizing  the  error  inherent  in  using 
a  time  delayed  DGPS  correction  was  developed.  It  consisted  of  simply  computing  the 
time  delay  differences  as  a  function  of  time  delay  and  fitting  those  differences  to  a  pre¬ 
sumed  model.  The  time  shift  ranged  from  zero  to  a  maximum  such  that  20  data  points 
remained  in  common  between  the  original  set  and  the  shifted  set.  Since  the  data  inter¬ 
val  was  10  seconds,  this  corresponds  to  an  overlap  of  200  seconds.  Thus  a  1500  sec- 
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ond  data  set  will  yield  a  1300  second  set  of  time  shifted  differences.  The  differences 
are  the  mean  square,  MS,  of  the  components  of  the  RMS  differences  at  each  time  shift. 

After  each  shift,  the  shifted  and  unshifted  errors  were  differenced.  The  differences 
(at  each  overlapping  point)  are  functions  of  the  total  shift,  T .  (A  shift  by  an  index  of 
one  corresponds  to  a  10  second  shift.)  For  the  discrete  case  these  differences  are 

=  Xerr.  ^-Xerr,  (3-9) 

1,  m  1  +  m  1 

The  mean  and  standard  deviation  at  each  shift,  m,  was  found  by  summing  over  i 


Miff;  . 


X.: 


—  _L 


(3.10) 


X  diff  =  — 


n„ 


(3.11) 


where  the  sum  on  i  is  over  the  number  of  overlapping  points,  n^,,  at  each  shift,  m. 
The  mean  square,  MS,  of  the  x,  y,  z  standard  deviations  was  computed 


dif4s(x)  =  x^diff^  +  y^diff„  +  z^diff„  (3.12) 

where  the  mean  square  difference  is  expressed  in  terms  of  the  time  delay,  t ,  corre¬ 
sponding  to  the  shift,  m. 
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Figure  3.3  shows  this  process  in  block  diagram  form.  The  result  of  the  time  delay 


error  process  is  a  curve  of  mean  square  error  for  differences  vs.  t  . 


Figure  3.3  Creating  Time  Delay  Errors 


3.6  Relation  of  Mean  Square  Differences  to  Markov  Model 

Consider  the  outputs  of  receiver  A  and  B:  (x^,  Ya,  Za)  and  (xg,  yg,  Zg).  Taking  the 

RSS  of  each  provides  RSSa  and  RSSg  as  in  Equation  (3.10).  The  difference  is  DIF'F  = 
RSSa  -  RSSg.  The  variance  on  the  difference  is  given  by 

®mFF  =  DIFF^-DIFF^  (3.13) 

If  the  difference  is  a  zero  mean  process,  then 

Oj-jjpp  =  DIFF  =  RSS^  —  2RSS^RSSg  +  RSSg  (3.14) 

Assume  both  A  and  B  are  zero  mean  processes,  then  RSS^  =  a^,  RSSg  =  Og,  and 
RSS^RSSg  =  po^Og .  This  leads  to 

^DIFF  " 


(3.15) 


2  2  2 

Assuming  A  and  B  have  similar  variances,  0^  =  0^  =  resulting  in 

Equation  (3.16). 


t  X  ^DIFF  “  "  ^^AB^^'P^ 

Letting  p  =  e  for  a  Markov  process  leads  to  the  model  function 


2 

^DIFF 
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1  -e 


t 

X 


X 

X 


c 


V  J 

which  is  equivalent  to  the  mean  square  error  if  the  bias  is  zero. 


(3.16) 


(3.17) 


3.7  Marquardt  Method 

The  Marquardt  Method  is  a  non-linear  curve  fitting  routine  that  enables  the  fit  of 
an  exponential  to  a  set  of  data.  This  is  the  means  of  obtaining  the  estimates  for  the 
DGPS  model  parameters.  The  inputs  to  the  model  include  arrays  of  the  independent 
and  dependent  data  points,  along  with  the  standard  deviations  for  each  data  point.  The 
algorithm  determines  the  standard  deviation  and  correlation  constant  for  the  first  order 
Markov  fit,  along  with  a  covariance  matrix  that  can  be  used  to  generate  error  bounds 
for  the  model. 

The  process  defines  a  chi-squared  merit  function  and  determines  best-fit 

parameters  by  minimizing  the  merit  function.  For  the  nonlinear  case,  the  minimization 

2 

is  performed  iteratively.  The  procedure  is  repeated  until  %  effectively  stops  decreas- 
ing.  The  selected  limit  for  the  difference  between  consecutive  %  values  was  0.1.  [22] 
The  Marquardt  method  is  described  in  more  detail  in  Appendix  A.  The  Marquardt 
method  is  used  to  determine  independently  the  errors  due  to  baseline  distance  and  time 
delay. 
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3.8  Pseudorange  Error  Processing 

All  previous  processing  was  performed  using  the  navigation  solutions  that  are  pro¬ 
vided  by  the  GPS  receiver.  Due  to  the  satellite  switches,  the  data  sets  are  limited  in 
length.  In  an  effort  to  obtain  longer  data  sets,  the  possibility  of  obtaining  pseudorange 
errors  was  investigated.  The  pseudorange  from  a  single  satellite  could  cross  over  sev¬ 
eral  satellite  switches  and  provide  a  source  of  continuous  data  for  use  in  the  process¬ 
ing. 

Equation  (3,18)  is  used  to  obtain  the  pseudorange  errors. 


P  -P 
-meas  -true 


=  H 


t  -r, 

meas  true 

At  -  At 

meas  true 


where  H  is  given  by  Equation  (3.19). 


(3.18) 


H  = 


T  1 

'I'  1 
»2  ^ 

1 

T 


(3.19) 


The  line  of  sight  vector  to  each  satellite  is  given  by  Ui^.  The  vector  p  contains 

-meas 

the  measured  pseudorange  for  each  of  the  four  satellites.  The  output  GPS  navigation 

solution  vector  is  r _ .  The  receiver’s  true  position  vector  is  r  .  The  user  clock 

bias  estimate  is  .  The  true  user  clock  bias  is  At  .  The  line  of  sight  vector  to 

each  satellite  is  obtained  by  using  post-processed  ephemeris  data  downloaded  from  the 
Internet.  The  ephemeris  data  is  obtained  at  15  minute  intervals,  and  an  interpolation 
routine  called  polint.m  is  used  to  find  the  true  ephemeris  for  all  data  points  [22]. 
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The  problem  with  the  pseudorange  processing  is  that  while  we  have  a  value  for 

rtrue»  the  true  clock  bias  in  Equation  (3.18)  is  not  known.  If  the  error  in  the  user  clock 

bias,  forced  to  be  zero,  the  pseudorange  errors  obtained  will  have 

jumps  in  them  at  each  satellite  switch.  This  is  shown  in  the  top  half  of  Figure  3.4.  If 

this  jump  is  attributed  to  the  change  in  At  due  to  the  satellite  switch,  a  correction 

mc8.s 

term  can  be  obtained  which  indeed  allows  the  three  continuing  pseudoranges  to  be 
continuous.  An  example  is  shown  in  the  bottom  half  of  Figure  3.4. 

The  actual  clock  error,  however,  is  still  unknown.  We  have  simply  found  the 
change  in  its  estimate  due  to  the  satellite  switch.  The  difference  between  the  measured 
and  true  clock  error  is  a  contributor  to  the  pseudorange  error,  and  cannot  be  ignored.  It 
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may  however  be  possible  to  estimate  the  contribution  of  this  term  statistically  based  on 
the  TDOP.  This  will  be  discussed  in  Section  5.2. 


3.9  Multipath 

As  GPS  signals  are  reflected  off  objects  before  entering  the  receiver  antenna,  a 
false  pseudorange  is  created.  This  happens  because  the  reflected  signal  has  traveled  a 
different  path  length  than  the  direct  signal.  As  noted  in  the  previous  section,  cyclic 
behavior  can  be  observed  in  the  GPS  position  output  which  may  be  attributable  to  this 
“multipath”  effect. 
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Oscillations  are  also  observed  in  the  time  delay  errors  formed  by  shifting  one 
receiver’s  data  in  time  and  differencing  this  with  the  other  receiver’s  output.  Figure  3.5 
illustrates  the  point  by  showing  the  time  delay  error  results  from  six  data  sets  in  the 
zero  baseline  case.  Note  the  oscillations  that  occur  in  each  set. 
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Figure  3.5  Oscillations  in  Time  Delay  Errors,  Zero  Baseline  Data 


After  noting  similar  frequencies  in  each  set,  the  average  period  was  identified  by 
taking  the  fast  Fourier  transform  (FFT)  of  each  set  shown.  Figure  3.6  shows  the  FFT 
outputs  from  each  of  the  six  sets.  The  FFT  output  is  very  similar  for  all  the  sets  with 
the  tallest  peaks  occurring  from  300  seconds  to  600  seconds.  Peaks  off  the  scale  were 
attributed  to  taking  the  FFT  of  a  finite  set  of  data.  Peaks  are  expected  to  occur  at  700  to 
1500  seconds  because  that’s  how  long  the  individual  data  segments  lasted.  The  pre- 
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dominant  periods  which  can  reasonably  be  attributed  to  multipath  are  between  300 
seconds  and  600  seconds. 
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Figure  3.6  FFT  of  Time  Delay  Errors,  Zero  Baseline  Data 


A  scenario  was  developed  which  results  in  the  multipath  frequencies  observed  for 


the  satellites  in  the  zero  baseline,  5-hour  data  set.  The  user  selects  a  position  offset 
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from  the  receiver  for  the  location  of  the  reflector.  The  true  range  and  satellite  position 
are  used  to  determine  the  multipath  frequency.  The  scenario  is  shown  in  Figure  3.7. 


Figure  3.7  Signal  Multipath 


As  shown,  dj  is  the  direct  path  from  satellite  to  receiver.  Two  possible  reflective 
paths,  d2,  were  considered  as  examples.  Note  that  d2  includes  the  path  before  and  after 
the  reflection  or  the  total  from  the  satellite  to  the  antenna.  One  was  a  reflection  off  a 
nearby  building.  The  position  offset  used  was  20  meters  north  and  20  meters  up.  The 
second  multipath  option  is  a  reflection  off  the  roof  that  enters  under  the  receiver.  This 
case  is  highly  likely  because  the  antenna  was  on  a  tripod  and  the  gain  under  the 
receiver  is  fairly  high.  The  offset  used  for  this  case  was  1.9  meters  down  and  0.5 
meters  north. 

The  multipath  frequency  is  given  in  Equation  (3.20)  and  Equation  (3.21). 


Ad,  Ad2 
A,At  A-At 


(3.20) 


^  d|(k+i)-d|(k)  d,(kti)-d,(k) 

0.19(t(k+ 1) -t(k))  0.19(t(k+l) -t(k))  '  ' 

where  X  is  the  wavelength  of  the  signal  (0.19  meters),  t  is  the  measurement  time,  and 

k  is  the  index  for  the  measurements.  The  results  for  satellite  26  for  the  first  offset  are 
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Multipath  Period  (seconds) 


shown  in  Figure  3.8.  The  period  is  between  90  and  100  seconds  (1.5  to  1.7  minutes) 
and  the  frequency  is  0.0095  to  0.011  Hz. 


GPS  Time  (hours) 

Figure  3.8  Sample  Multipath  Period  for  Offset  of  20  Meters 


Figure  3.9  shows  the  results  for  the  reflections  under  the  antenna.  The  period  in  this 
case  is  150  seconds  or  greater  (1.5  - 10  minutes),  and  the  frequency  is  0.008  Hz  or  less. 
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Figure  3.9  Sample  Multipath  Period  for  Offset  of  1.9  Meters 


Experimental  observations  of  multipath  were  found  in  [1,  16].  Multipath  frequen¬ 
cies  with  time  periods  of  5  to  30  minutes  were  recorded.  Differential  multipath  with 
time  periods  of  5  to  10  seconds  were  observed.  The  simulation  results  are  in  the  range 
of  these  observations.  Other  multipath  “paths”  are  possible.  The  two  examples  were 
only  intended  to  indicate  reasonableness.  No  attempt  was  made  to  actually  measure 
the  distance  to  nearby  buildings. 

Note  that  the  frequencies  observed  in  the  data  and  those  found  from  the  simulation 
of  reflections  from  under  the  receiver  fall  into  the  same  frequency  range.  The  fre¬ 
quency  for  reflections  off  nearby  buildings  matches  closer  with  the  frequencies  found 
in  the  raw  data  highlighted  in  the  last  section. 
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Chapter  4 


Results 

4.1  Rooftop  Experiment  Results 

4.1.1  Navigation  Solution  for  5-Hour  Data 

The  first  set  of  results  is  from  processing  the  navigation  solution  for  the  5-hour  set 

of  data  with  the  zero  baseline.  Figure  4.1  shows  the  satellites  that  were  used  during  the 
data  collection.  The  y-axis  is  the  satellite  identification  number.  Note  that  at  each  time 
only  four  satellites  are  being  used.  The  gap  between  115  and  116.5  hours  was  noted 
earlier  as  the  time  when  less  than  four  satellites  were  tracked. 


Figure  4.1  Satellite  ID  vs.  Time 


The  vertical  lines  mark  each  spot  where  there  is  a  switch  in  at  least  one  of  the  sat¬ 
ellites  being  tracked.  The  data  was  split  into  six  different  sets  where  the  satellite  set 
stayed  constant  over  that  interval.  The  small  data  segment  towards  the  middle  was  not 
used. 
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Time  delay  differences  were  created  using  the  processing  described  in  Chapter  3. 
Figure  4.2  shows  the  resulting  mean  square  vs.  time  delay  plots  for  each  of  the  six  sets 
of  data.  The  average  is  also  plotted  on  this  chart. 


Figure  4.2  Navigation  Solution  Errors 


Figure  4.3  shows  the  1-sigma  bounds  on  the  time  delay  differences.  These  were 
computed  by  finding  the  standard  deviation  for  the  mean  square  error  calculation  as 
shown  in  Equation  (4.1).  [11]  Note  that  the  error  bounds  for  the  six  sets  overlap  much 
of  the  time. 
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2 

where  is  the  mean  square  error  (ordinate  in  graph),  and  n^,  is  the  number  of 
overlapping  points.  It  can  be  seen  that  the  error  bounds  increase  as  the  number  of  over¬ 
lapping  points  decreases  for  increasing  time  delay. 


Figure  4.4  shows  the  weighted  average  from  the  six  data  sets,  along  with  the  stan¬ 
dard  deviation  on  the  average.  These  values  are  used  as  inputs  to  the  Marquardt  curve 
fitting  routine.  The  weighted  average  and  standard  deviation  were  found  using  Equa¬ 
tion  (4.2)  and  (4.3)  where  MSj  is  the  mean  square  for  each  set  and  Gj  is  the  standard 
deviation  for  the  mean  square. 


mean 


(4.2) 
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Figure  4.5  shows  the  Markov  fit  to  the  average  from  the  6  data  sets.  Equation  (4.4) 
shows  the  model  that  was  used  in  the  Marquardt  fitting  routine.  A  sum  of  two  Markov 
processes  was  used  because  this  model  provided  a  better  fit  to  the  data.  The  first  pro¬ 
cess  with  the  shorter  time  constant  is  attributed  to  multipath.  The  period  falls  in  the 
sample  range  found  in  the  multipath  scenario  as  discussed  in  Chapter  3.  All  the  zero 
baseline  cases  exhibit  this  short  time  constant  behavior.  The  longer  time  constant  pro¬ 
cess  is  intended  to  encompass  all  other  errors.  A  model  for  this  process  was  one  of  the 
goals  of  the  study. 


Model  = 
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+  02 
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The  error  bounds  on  the  Markov  fit  are  shown  with  dotted  lines.  These  are 


obtained  by  using  the  covariance  provided  by  the  Marquardt  routine  as  shown  in 
Equation  (4.5). 


s  afcSf 


(4.5) 


The  one-sigma  error  on  the  model  is  8e .  The  function  f  is  the  mean  square  error 
which  is  a  function  of  the  vector  s .  The  vector  s  consists  of  both  values  for  a  and  i . 
The  partial  derivatives  of  the  function  f  with  respect  to  the  elements  of  s  are  output 
from  the  Marquardt  routine.  The  routine  also  provides  the  covariance  matrix  E  defined 
by  Equation  (4.6). 


T 

E  =  8s5s  (4.6) 

The  noise  level  for  the  difference  between  receiver  A  and  B  is  0.05  meters^  and 
provides  a  bias  for  the  navigation  solution  errors.  The  receiver  noise  has  a  mean  of 
0.21  meters  and  a  standard  deviation  of  0.23  meters.  Squaring  the  standard  deviation 
gives  a  bias  for  the  mean  square  of  0.05  meters^. 

Figure  4.6  shows  the  residual  for  the  Markov  process  fit.  This  is  obtained  by  sub¬ 
tracting  the  fit  from  the  data.  The  residual  represents  anything  that  has  not  been 
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Residual 


bias  =  0.05  meters 


included  in  the  model.  Note  the  sinusoidal  characteristic  in  the  residual.  This  will  be 


addressed  later  in  the  chapter. 


Figure  4.6  Residual  for  Markov  Fit 
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The  model  parameters  for  each  of  the  6  sets  and  the  average  are  given  in  Table  4.1. 
The  average  values  are  obtained  from  a  separate  fit  for  the  average  curve,  not  just  the 
average  for  the  parameter  values  from  the  6  data  sets.  The  fits  to  the  individual  sets 
reveal  the  fact  that  the  Markov  model  is  not  very  sensitive  to  changes  in  the  X2  value. 
That  is  the  reason  for  the  large  variations  seen  in  this  parameter.  This  insensitivity  is 
related  to  the  limit  (variance)  of  the  experimental  correlation  given  in  Chapter  3. 


Table  4.1  Results  for  Zero  Baseline 


SET 

(meters) 

'^1 

(seconds) 

<^2 

(meters) 

'^2 

(seconds) 

1 

1.25 

272.5 

285.9 

2 

0.68 

80.2 

5.74 

7643.3 

3 

0.85 

39.8 

26.4 

123,719 

4 

1.86 

204.9 

5.80 

11,412 

5 

1.47 

125.2 

3.92 

10,110 

6 

1.10 

57.6 

506.6 

Average 

0.98 

103.9 

1.71 

499.0 

4.1.2  Navigation  Solution  for  24-Hour  Data 

The  second  set  of  results  is  from  processing  the  navigation  solution  for  the  24-hour 

set  of  data  with  the  zero  baseline.  Figure  4.7  shows  the  satellites  that  were  used  during 
the  data  collection.  The  y-axis  is  the  satellite  identification  number.  Note  that  at  each 
time  only  four  satellites  are  being  used.  The  vertical  lines  mark  each  spot  where  there 
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ID  Number 


is  a  switch  in  at  least  one  of  the  satellites  being  tracked.  Ten  sections  were  selected  for 
processing  of  the  navigation  solution.  These  are  indicated  by  an  ‘X’  above  the  data  set. 


Time  delay  differences  were  created  and  Figure  4.8  shows  the  resulting  mean 
square  vs.  time  delay  plots  for  each  of  the  ten  sets  of  data.  The  average  is  also  plotted 
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on  this  chart.  The  average  is  the  very  thin  line  that  runs  through  the  middle  of  the 


curves. 


Noting  that  some  of  the  curves  in  Figure  4.8  are  lower  than  others,  the  time  of  day 
for  each  curve  and  the  mean  MSE  were  calculated.  Figure  4.9  plots  the  mean  of  each 
of  these  sets  of  time  delay  differences  versus  time  of  day.  The  time  of  day  is  plotted 
from  0  to  24  hours  (1  days).  Hour  zero  occurs  at  midnight.  It  might  be  said  that  there  is 
a  large  mean  error  between  8:00  a.m.  and  10:00  a.m.  corresponding  to  a  rise  in  total 
electron  count  in  the  morning,  the  smallest  mean  square  differences  occur  between 
noon  and  8:00  p.m.  There  is  then  another  increase.  This  increase  occurs  significantly 
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after  sunset.  Any  correlation  between  the  mean  square  differences  and  sunrise  and 
sunset  is  difficult  to  discern. 


Time  of  Day  (Hours) 

Figure  4.9  Mean  Error  by  Time  of  Day 


Figure  4.10  shows  the  average  from  the  ten  data  sets,  along  with  the  standard  devi¬ 
ation  on  the  average.  These  values  are  used  as  inputs  to  the  Marquardt  curve  fitting 
routine. 


Figure  4.10  Average  Navigation  Solution  Error,  24-Hour  Data 
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Figure  4.11  shows  the  Markov  fit  to  the  average  from  the  10  data  sets.  The  bias  is 
zero  because  there  is  only  one  receiver  for  the  24-hour  data.  The  result  is  a  shift  of 
receiver  A  vs.  receiver  A.  The  T2  value  of  3847.1  seconds  equals  1.07  hours. 


Note  that  the  error  bounds  on  the  raw  data,  Figure  4. 10,  and  on  the  Markov  fit,  Fig¬ 
ure  4.11,  do  not  completely  overlap.  This  is  especially  obvious  after  a  few  thousand 
seconds.  This  “inconsistency”  reflects  the  fact  that  the  Markov  process  is  an  incom¬ 
plete  description  of  the  pseudorange  transmission  delay  differences.  In  fact,  the  errors 
and  their  differences  are  not  stochastic,  but  are  “unpredictable”.  Even  this  is  not 
exactly  the  case  because,  with  enough  effort,  satellite  ephemeris  and  transmission 
delay  could  be  predicted  at  least  in  the  short  term.  Characterizing  these  differences  as 
Markov  processes  is  merely  convenient.  As  the  analysis  proceeds  the  character  of  the 
fit  residuals  will  be  observed. 
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Figure  4.12  shows  the  residual  for  the  Markov  fit  on  the  24-hour  data.  Again  note 
the  sinusoidal  nature  of  the  model  error. 


Figure  4.13  shows  the  Markov  fit  for  both  the  5-hour  and  24-hour  data  sets.  The 
model  parameters  for  these  two  sets  are  almost  identical.  The  data  was  collected  on 
different  days  and  times  of  day.  The  results  indicate  a  definite  consistency  between  the 
data  sets.  Note  that  at  the  end  of  each  data  set  the  average  error  increases  sharply.  This 
is  a  result  of  the  fact  that  the  number  of  data  points  is  decreasing  at  each  shift.  The 
smaller  segments  have  only  20  points  left  at  the  end,  or  only  200  seconds  of  data.  It 
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Mean  Square  Error  (melers''2) 


seems  reasonable  that  there  would  be  trouble  estimating  a  1-hour  time  constant,  with 
so  few  points  remaining. 


Figure  4.13  Markov  Process  Fit,  Both  5-Hour  and  24-Hour  Data 


The  fit  to  each  of  the  10  navigation  solution  sets  is  given  in  Table  4.2.  The  large 
variance  on  the  Xj  values  is  again  apparent,  and  we  will  note  once  more  that  it  is  con¬ 
sistent  with  the  variance  in  the  correlation  given  in  Chapter  3. 


Table  4.2  Results  for  Zero  Baseline,  24-Hour  Data 


SET 

(meters) 

(seconds) 

(meters) 

^2 

(seconds) 

1 

1.57 

89.8 

9.51 

38,358 

2 

1.81 

255.0 

12.54 

45,316 

3 

1.50 

208.6 

1.27 

214.0 

4 

1.07 

123.1 

41.4 

7.4x10^ 

5 

0.93 

105.1 

15.3 

46,620 

6 

1.98 

154.3 

1.36 

512.5 

7 

2.41 

296.7 

3.76 

22,579 

8 

0.69 

43.4 

13.1 

67,439 

9 

1.77 

252.0 

0.90 

261.6 
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Table  4.2  Results  for  Zero  Baseline,  24-Hour  Data 


SET 

(meters) 

(seconds) 

a2 

(meters) 

(seconds) 

10 

1.54 

120.1 

32.1 

2.7x10^ 

Average 

1.47 

171.8 

2.77 

3847.1 

4.1.3  Pseudorange  Errors,  5-Hour  Data 

The  third  set  of  results  is  from  processing  the  pseudorange  errors  for  the  5-hour  set 

of  data  with  the  zero  baseline.  Fourteen  pseudorange  sequences  were  formed  from  the 
data.  The  large  gap  in  the  data  forced  the  splitting  of  pseudoranges  on  either  side  of  the 
gap.  Figure  4.14  shows  the  resulting  mean  square  vs.  time  delay  plots  for  each  of  the 
fourteen  sets  of  data.  The  average  is  also  plotted  on  this  chart. 


Pseudorange  Error  vs  Time  Delay,  O  Baseline 
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Mean  Square  Error 
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Residual  (meter''2) 


The  residual  for  the  Markov  fit  is  shown  in  Figure  4.17.  Once  again  note  the  sinu¬ 
soidal  nature  of  the  model  error. 


Figure  4.17  Residual  for  Markov  Fit,  Pseudorange  Errors 


4.1.4  Pseudorange  Errors,  24-Hour  Data 

The  fourth  set  of  results  is  from  processing  the  pseudorange  errors  for  the  24-hour 
set  of  data  with  the  zero  baseline.  Ten  pseudorange  sequences  were  formed  from  the 


68 


data.  Figure  4.18  shows  the  resulting  mean  square  vs.  time  delay  plots  for  each  of  the 
ten  sets  of  data.  The  average  is  also  plotted  on  this  chart. 
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Time  Delay  (seconds) 

Figure  4.18  Pseudorange  Error,  24-Hour  Data 


Figure  4.19  shows  the  average  from  the  fourteen  data  sets,  along  with  the  standard 
deviation  on  the  average.  These  values  are  used  as  inputs  to  the  Marquardt  curve  fit¬ 
ting  routine. 
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Figure  4.20  shows  the  Markov  fit  to  the  average  from  the  14  data  sets.  Again,  The 


first  rise  is  very  close  to  what  was  seen  in  the  navigation  solution  results. 


Time  Delay  (seconds) 


Figure  4.20  Markov  Fit  for  Pseudorange  Error 


The  residual  for  the  Markov  fit  is  shown  in  Figure  4.21.  Once  again  note  the  sinu¬ 
soidal  nature  of  the  model  error. 


Time  Delay  (seconds) 

Figure  4.21  Residual  for  Markov  Fit,  Pseudorange  Errors 
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Mean  Square  Error  (meters*2) 


The  Markov  fit  for  all  cases  is  shown  in  Figure  4.22.  All  the  curves  are  aligned  in 
the  first  rise. 


4.1.5  Residual  Errors 

All  the  residual  errors  shown  have  had  sinusoidal  characteristics.  They  are  now 
examined  more  closely.  Figure  4.23  plots  the  model  residuals  for  all  three  of  the  previ- 
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ous  cases.  Not  only  are  the  sinusoidal  characteristics  evident  in  each,  but  the  same  fre¬ 
quency  and  phase  is  present. 

0.5 
0.4 
0.3 
0.2 
0.1 
0 

-0.1 
-0.2 
-0.3 
-0.4 
-0.5 

O  100  200  300  400  500  600  700  800 

Time  Delay  (seconds) 

Figure  4.23  Residual  Errors,  All  Cases 

Figure  4.24  shows  the  FFT  for  each  of  the  model  residuals.  The  periods  of  the 
oscillations  are  very  similar  for  all  four  cases.  A  primary  period  of  about  400  seconds 
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is  apparent.  The  next  step  is  to  see  if  the  Markov  model  fits  the  data  better  if  appropri¬ 
ate  sinusoidal  terms  are  added. 


o  100  200  300  400  500  600 


Period  (sec) 


Figure  4.24  FFT  for  Residuals 

A  sinusoidal  term  using  the  period  identified  was  added  to  the  Markov  model  to 
attempt  a  better  fit.  The  new  fit  is  shown  in  Figure  4.25.  This  fit  is  much  better  than  the 
Markov  process  alone.  The  model  for  the  residual  used  is  given  in  Equation  (4.7).  The 
short  period  of  this  oscillation  once  again  points  to  multipath.  Some  particular,  strong 
reflected  signal  common  to  all  cases  is  a  possible  explanation 


Mean  Square  Error  (meters*2) 


Res  =  0.35sin  (t  +  150) 


(4.7) 


Figure  4.25  New  Model  Fit  to  Data 


Figure  4.26  shows  the  residual  for  the  new  model.  Note  that  the  amplitude  of  the 
residual  has  been  reduced  to  about  a  fourth  of  its  former  level.  Note  that  sinusoidal 
characteristics  can  still  be  observed  in  the  data.  Further  reduction  using  these  periods 
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is  possible.  The  primary  periods  in  these  residuals  match  those  found  in  the  time  delay 
data  in  Chapter  3. 


;  East  and  West  Coast  Experiment  Results 

Figure  4.27  shows  the  results  for  the  East  Coast  baselines. 


Figure  4.28  shows  the  results  from  the  West  Coast  data. 


The  errors  shown  in  Figure  4.27  and  Figure  4.28  include  the  effects  of  multipath. 
An  estimate  of  errors  due  to  multipath  was  obtained  using  the  same  process  that  was 
used  for  the  zero  baseline  data.  That  is,  time  shift  differences  were  made,  and  a  dual 
Markov  model  was  fit  to  these  differences.  The  shorter  time  constant  was  consistent 
with  multipath  effects.  The  variances  were  somewhat  smaller  than  seen  on  the  rooftop. 
Variances  for  each  pair  were  obtained  by  combining  the  appropriate  values.  These 
variances  are  shown  in  Table  4.3.  They  were  subtracted  from  the  errors  before  pro¬ 
ceeding  to  find  a  combined  time  delay/distance  model. 


Table  4.3  Multipath  Variances  for  Each  Baseline 


Separation  (km) 

Locations 

(meters^) 

84 

Draper/Portsmouth 

0.22 

149 

Long  Beach/San  Diego 

0.92 

282 

Fallon/Stockton 

0.67 

496 

Fallon/Fort  Irwin 

0.90 
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Mean  Square  Error  (meters''2) 


Table  4.3  Multipath  Variances  for  Each  Baseline 


Separation  (km) 

Locations 

(meters^) 

740 

Langley/Draper 

0.18 

756 

Fallon/San  Diego 

0.49 

820 

Langley/Portsmouth 

0.10 

4.3  Combined  Model 

The  results  have  been  presented  for  each  of  the  three  data  sources:  zero  baseline, 
East  Coast,  and  West  Coast.  Now  the  results  from  the  navigation  solution  errors  will 
be  combined  to  create  composite  results.  The  24-hour  data  for  the  zero  baseline  case 
will  be  used  since  it  has  a  longer  time  base. 

4.3.1  Results  From  All  Sets 

Figure  4.29  plots  all  baseline  errors  together. 


o  200  400  600  800  1 000  1 200 

Time  Delay  (seconds) 


Figure  4.29  Navigation  Solution  Errors,  All  Data  Sets 


Figure  4.30  plots  the  same  information  in  3-D  format  for  a  more  visual  interpreta¬ 


tion. 


Baseline  Distance  (km)  0  0  (seconds) 

Figure  4.30  3-D  Error  Plot,  All  Data  Sets 


Figure  4.31  plots  a  cross  section  of  the  3-D  plot  at  zero  time  delay.  This  cross  sec¬ 
tion  will  be  used  to  fit  the  Markov  process  for  the  distance  correlation.  The  last  three 
points  are  from  the  East  Coast  baselines  with  a  C/A-code  receiver  at  one  end  (Willow 
Run).  The  Draper  data  collection  personnel  indicated  a  lower  confidence  in  the  output 
from  these  receivers.  In  addition,  the  selective  aveiilability  errors  had  to  removed  by 
Aerospace  in  post  processing.  By  policy,  these  errors  were  not  totally  removed  from 
the  data,  and  some  noise  still  remains.  Given  these  facts,  a  strong  case  can  be  made  for 


78 


Mean  Square  Error  (meters''2) 


4.3.2  Results  From  High  Confidence  Sets 

To  create  a  set  of  high  confidence  data,  the  three  C/A-code  sets  at  838,  1022  and 

1052  km  were  removed.  The  resulting  time  delay  errors  are  shown  in  Figure  4.32. 


Figure  4.32  Time  Delay  Errors,  High  Confidence  Sets 
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Figure  4.33  shows  the  3-D  plot  for  the  high  confidence  data  sets.  Note  that  this  sur¬ 


face  is  smoother  than  the  one  including  the  C/A-code  results. 


Baseline  Distance  (km)  0  0  (seconds) 

Figure  4.33  3-D  Error  Plot,  High  Confidence  Sets 


Figure  4.34  shows  the  3-D  plot  for  the  high  confidence  data  sets  with  the  multipath 
removed  from  the  zero  baseline  data.  This  plot  will  be  compared  with  the  final  model. 


Figure  4.34  3-D  Error  Plot,  No  Multipath 


81 


Figure  4.35  plots  the  cross-sectional  view  of  Figure  4.33  for  the  case  with  no  time 
delay.  This  data  is  used  as  an  input  to  the  Marquardt  method  to  fit  the  Markov  model  to 
the  distance  dimension.  The  1 -sigma  error  bounds  provided  to  the  Marquardt  method 
are  also  shown  in  the  graph.  The  circles  show  the  actual  data  points  available  for  the 
Marquardt  fit.  Note  that  even  with  several  baselines  of  collected  data,  there  are  still 
very  few  points  in  the  distance  dimension  compared  with  the  number  of  points  in  the 
time  delay  dimension. 


Figure  4.36  shows  the  Markov  fit  for  the  distance  dimension.  The  bias  is  zero 
because  the  24-hour  data  was  used  at  0  km,  and  only  one  receiver  was  available  for 
this  data.  Equation  (4.8)  shows  the  model  that  was  used  in  the  Marquardt  fitting  rou¬ 
tine  for  the  distance  dimension. 
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Mean  Square  Error  (meters^2) 


Baseline  Distance  (km) 

Figure  4.36  Markov  Process  Fit,  Distance  Dimension 


Figure  4.37  plots  the  residual  for  the  distance  dimension  Markov  fit.  There  is  a 
slight  sinusoidal  characteristic  to  this  plot,  but  not  enough  data  points  are  available  to 
merit  trying  to  correct  the  model. 


Figure  4.38  shows  the  final  resulting  DGPS  model  obtained  by  combining  the 
Markov  fit  to  the  time  delay  and  baseline  distance  dimensions.  Since  the  same  under¬ 
laying  phenomena  cause  the  distance  and  time  decorrelation,  the  variance  of  the  pro¬ 
cess  is  common.  The  value  of  3.73  meters^  is  the  weighted  average  of  the  values  found 
for  the  time  and  distance  models.  Note  that  the  raw  data  has  a  significant  spread  and 
this  model  will  not  encompass  every  case.  Rather  it  is  meant  to  provide  an  average 
error  model  that  represents  the  typical  behavior  of  the  DGPS  operation.  The  formula 
for  the  DGPS  model  is  given  by  Equation  (4.9).  This  plot  can  be  compared  with  Figure 
4.34.  Note  that  these  results  have  been  normalized  by  the  DOP  values.  To  determine 
the  errors  obtained  in  the  navigation  solution,  the  model  must  be  multiplied  by  the 
appropriate  DOP  values  for  the  given  satellite  geometry. 
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Figure  4.38  DGPS  Error  Model 

The  data  used  to  determine  the  time  constant,  x ,  comes  exclusively  from  the  ‘‘roof¬ 
top”  experiments.  This  is  because  the  data  from  the  East  and  West  Coast  experiments 
was  not  of  sufficient  duration  to  contribute  much  to  the  determination  of  T . 
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Chapter  5 


Conclusions  and  Suggestions  for  Future  Work 


This  thesis  analyzed  GPS  data  collected  by  Draper  Laboratory  to  determine  DGPS 
error  characteristics.  It  was  found  that  a  first  order  Markov  process  could  be  used  as  an 
appropriate  model  for  the  errors  in  differential  GPS.  Since  the  errors  are  correlated  in 
both  time  and  distance,  the  Markov  model  has  a  process  variance,  a  time  correlation 
constant,  and  a  distance  correlation  constant.  Analysis  of  both  navigation  solution 
errors  and  pseudorange  bias  errors  provided  a  wide  range  of  data  sets  for  processing. 
The  Marquardt  non-linear  curve  fitting  routine  was  used  to  obtain  the  parameters  for 
the  Markov  model  for  each  set  of  data. 


5.1  Conclusions 

The  resulting  DGPS  model  is  given  in  Equation  (5.1). 


MS  =  3.73 


1-e 


t _ X  \ 

3847.1  122.8 


meters 


(5.1) 


V  J 

where  MS  is  the  mean  square  error  in  the  pseudorange  biases. 

This  model  is  fairly  representative  of  average  navigation  solution  errors.  It  should 
be  noted  that  there  are  significant  deviations  from  nominal  performance.  (Refer  to  Fig¬ 
ure  4.8  for  example.)  The  analyst  might  wish  to  choose  a  larger  value  for  the  variance 
to  describe  “worst  case”  performance.  The  rationale  for  using  a  Markov  fit  is  that 
many  physical  processes  tend  to  exhibit  Markov-like  behavior  and  that  it  is  a  conve¬ 
nient  model  for  linearized  analysis. 
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The  error  model  is  consistent  with  the  differential  GPS  error  budget  given  in  Chap¬ 
ter  1.  The  DGPS  budget  had  a  MS  error  of  1.59  meters^,  (1.26^).  The  Markov  model 
determined  herein  has  a  MS  error  of  3.73  meters^.  The  variance  at  92.6  km  separation 
and  zero  time  delay  is  about  1.98  meters^.  Fifty  nautical  miles  (92.6  kilometers)  is  the 
distance  assumed  in  the  referenced  error  budget.  Multipath  effects  could  easily 
account  for  the  difference.  The  referenced  budget  is  only  the  specification  for  the  sys¬ 
tem  and  not  necessarily  indicative  of  actual  performance. 

The  raw  data  has  a  significant  spread  to  it  depending  on  the  time  of  day,  geometry, 
multipath,  weather  conditions,  and  possibly  other  factors.  This  model  does  not  encom¬ 
pass  every  facet  of  the  error  sources.  In  particular  we  have  seen  the  importance  of  mul¬ 
tipath  error  modeling.  An  attempt  was  also  made  to  correlate  the  magnitude  of  the 
time  delay  errors  to  the  time  of  day.  Few  data  points  were  available,  but  the  results 
indicate  the  lowest  values  occur  during  the  afternoon  and  early  evening.  This  would 
mean  that  the  output  from  the  two  receivers  is  more  highly  correlated  at  these  times. 

The  model  was  developed  by  analyzing  the  errors  in  the  navigation  solution  output 
of  the  GPS  receivers.  The  individual  pseudorange  errors  that  contribute  to  the  position 
error  were  also  examined.  A  technique  was  developed  to  attribute  changes  in  all  pseu¬ 
dorange  errors  at  satellite  switch  to  changes  in  the  user  clock  bias  error  estimate.  This 
indeed  caused  the  pseudorange  error  for  each  satellite  to  be  continuous  and  not  experi¬ 
ence  jumps  when  a  satellite  switch  occurred.  The  results  from  this  analysis  were  con¬ 
sistent  with  those  obtained  from  the  navigation  solution  analysis  for  the  zero  baseline 
data.  While  this  is  encouraging,  and  lends  more  credibility  to  the  determination  of  the 
time  constant,  the  lack  of  a  measure  for  true  GPS  time  means  the  value  for  pseudor¬ 
ange  error  determined  in  this  manner  will  be  consistently  low. 
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The  cross-correlation  technique  was  ineffective  at  determining  the  Markov  fit 
parameters.  This  is  due  to  the  requirement  for  large  continuous  data  sets  for  cross-cor¬ 
relation  accuracy.  DGPS  data  sets  are  limited  to  about  one  hour  due  to  satellite 
switches.  Alternatively,  the  Marquardt  curve  fitting  routine  was  used  to  fit  the  time 
delay  error  to  a  Markov  process.  This  routine  also  allowed  fitting  of  two  independent 
processes. 

This  thesis  effort  has  revealed  many  interesting  aspects  of  the  GPS  receiver  opera¬ 
tion  including  a  close  look  at  P-code  data  from  two  receivers  attached  to  the  same 
antenna.  Of  biggest  importance  may  be  the  dramatic  change  in  amplitude  and  fre¬ 
quency  of  oscillations  in  the  navigation  solution  at  each  point  where  a  satellite  switch 
occurs.  Given  that  multipath  is  the  cause  of  such  behavior,  shielding  or  other  counter¬ 
measures  should  be  investigated  on  future  experiments  or  in  operational  scenarios. 

5.2  Future  Work 

As  just  mentioned,  multipath  should  be  treated  more  carefully  in  future  experi¬ 
ments.  It  should  be  eliminated  in  order  to  examine  other  error  sources  and  should  be 
introduced  in  a  controlled  fashion  if  desired. 

Internet  GPS  data  may  be  able  to  provide  longer  sets  of  data  where  the  restriction 
of  common  satellites  in  each  set  is  no  longer  a  requirement.  In  order  to  obtain  this,  the 
best  sites  with  accurately  surveyed  positions  and  hydrogen  maser  clocks  would  pro¬ 
vide  the  best  opportunity  for  such  data.  This  data  would  provide  a  check  against  the 
work  done  in  this  thesis  and  provide  more  insight  into  the  correlation  constants 
obtained.  The  best  Internet  address  to  use  for  these  sources  is  toba.ucsd.edu.  It  is  not 
clear,  however,  how  the  errors  such  as  selective  availability  can  be  completely  sepa¬ 
rated  from  the  other  physical  errors  of  interest.  The  Internet  receivers  are  not 


89 


encrypted,  so  SA  will  be  present.  Since  the  position  is  known  very  accurately,  the 
errors  in  the  navigation  solution  will  be  known  very  accurately,  but  the  contribution 
from  each  error  source  can  only  be  estimated. 

Even  if  true  GPS  time  is  unavailable  at  the  receiver,  a  comparison  of  Time  Dilution 
of  Precision  (TDOP)  values  with  the  PDOP  values  could  be  utilized  to  obtain  a  statis¬ 
tical  estimate  of  the  receiver  clock  errors.  This  would  allow  inclusion  of  these  errors 
into  the  estimate  of  the  pseudorange  errors  combined  with  the  adjustment  of  clock  bias 
devised  herein.  This  would  allow  longer  data  sets  to  be  analyzed.  Short  segments  of 
data  were  one  of  the  primary  limitations  for  this  thesis  effort. 
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Appendix  A 


Marquardt  Method 

The  Marquardt  method  or  Levenberg-Marquardt  method  provides  a  way  to  per¬ 
form  non-iinear  curve  fitting  for  data.  The  process  defines  a  chi-squared  merit 
function  and  determines  best-fit  parameters  by  minimizing  the  merit  function.  For  the 

nonlinear  case,  the  minimization  is  performed  iteratively.  The  procedure  is  repeated 

2 

until  X  effectively  stops  decreasing.  Numerical  Recipes  [22]  explains  the  Marquardt 
method  in  detail.  The  main  points  will  be  discussed  here. 

The  model  to  be  fitted  is 


y  =  y  (x^a) 


(A.2) 


where  a  is  the  vector  of  parameters  in  the  model  that  the  Marquardt  method  must  find 

r  xA 


In  this  case,  a  =  |a  and  y  =  & 


1  -e 

V  J 


.  The  X  merit  function  is 


N  ^ 


i=  1 


Yi-y  (Xi;a) 


<7; 


(A.3) 


The  goal  of  the  Marquardt  method  is  to  find  a^in  that  minimizes  x  •  An  initial  guess 
called  a  cur  is  required  to  get  things  started.  If  the  guess  for  a  is  fairly  good,  the  second 
derivative  matrix  (Hessian  matrix)  of  the  x  function  is  used  to  jump  immediately  to 
the  a  for  that  guess  using 


^min  =  +  *  ["Vx^Ca^j] 


(A.4) 


where  D  is  the  Hessian  matrix. 
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However,  if  the  guess  is  poor,  the  inverse-Hessian  method  cannot  be  used.  In  this 
case,  the  steepest  descent  method  is  used  to  take  a  step  down  the  gradient 

^next  =  ^cur "  Constant  x  (a^„^)  (A.5) 

where  the  constant  is  small  enough  not  to  exhaust  the  downhill  direction. 

2  _  ■  2 
The  gradient  of  %  with  respect  to  the  parameters  a ,  which  will  be  zero  at  the  % 

minimum,  has  components 


Taking  an  additional  partial  derivative  gives 


(Xi;a) 


k  =  1,2,...,M 


(A.6) 


3^  = 


(A.7) 


The  following  terms  are  defined 


3 

2aa. 


Km  =  X- 


2  9 

_  z 


‘^‘■2aa,aa, 


(A.8) 


Now  [a]  =  -D  in  Equation  (A.4)  which  can  be  rewritten 


M 

^a^i5a,  =  3,^  (A.9) 

i=  1 

This  set  of  linear  equations  can  be  solved  for  the  increments  5  a,  that  are  added  to  the 
current  approximation  to  get  the  next  approximation. 

The  steepest  descent  formula  translates  to 


5a,  =  constant  X  3, 


(A.  10) 
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The  Marquardt  method  allows  the  smooth  variation  between  the  steepest  descent 
method  and  the  inverse-Hessian  method  as  the  approximation  goes  from  poor  to  good. 
The  two  methods  are  combined  into  one  equation 

M 

^aki6a,  =  Pi,  (A.11) 

i=  1 


where  a  is  given  by 


ay^a^-H+X) 


(A.12) 


and  X  is  referred  to  by  the  author  as  a  non-dimensional  fudge  factor.  If  X  is  large,  the 
approximation  is  poor,  and  the  steepest  descent  method  is  used.  If  X  is  near  zero,  the 
approximation  is  good,  and  the  inverse-Hessian  method  is  used. 


The  Marquardt  recipe  consists  of  the  following: 

2  _ 

•  Compute  X  (a) 

•  Pick  a  modest  value  for  X.,  say  =  0.001. 

•  (1)  Solve  Equation  (A.  11)  for  5a  and  evaluate  %  (a  -i-  5a) 

2  2 

•  If  X  (a  +  5a)  >  X  (a)  ,  increase  X  by  a  factor  of  10  and  return  to  (1) 

2  2 

•IfX  (a  +  5a)<x  (a),  decrease  X  by  a  factor  of  10,  update  the  trial  solution 
a  a  +  5a ,  and  return  to  (1). 

2 

Once  the  acceptable  minimum  is  found  (change  in  %  of  less  than  0.1),  X  is  set  to 
zero  and  the  covariance  is  found 


[C]  =  [a]“‘  (A.13) 

The  partial  derivatives  of  the  model  with  respect  to  the  individual  parameters  in  a 
are  used  by  the  Marquardt  method.  For  the  Markov  exponential  fit,  the  partials  are 


^  = 
d<5 


2c 


X 

) 


(A.14) 
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(A.  15) 


9x 


_x 

2  T 

-a  xe 

2 

T 
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Appendix  B 


Autocorrelation  Limits 


The  limits  on  the  experimental  autocorrelation  function  were  found  in  [7].  The  ini¬ 
tial  assumption  proven  in  [3]  is  that  the  variance  of  an  experimentally  determined 
autocorrelation  function  satisfies  the  inequality 


VarV^(T)  <|jR^(i:)dx 


where  T  =  time  length  of  the  experimental  record 

(x)  =  autocorrelation  function  of  the  Gaussian  process  under  consideration 
Vjj  (x)  =  time  average  of  X-r  (t)  Xy  (t  +  x)  where  Xy(t)  is  the  finite-length  sam¬ 
ple  of  X  (t)  (that  is,  (x)  is  the  experimentally  determined  autocorrelation  function 
based  on  a  finite  record  length) 

The  experimental  autocorrelation  is  formed  by 


Vx('c)  =  ^J^^"^XT-(t)XT(t  +  x)^ 


It  is  noted  here  that  the  Matlab  function  xcorr  does  not  divide  by  =7^ ,  so  the  user 

T-x 

must  be  careful  when  working  with  that  function.  Dividing  properly  will  result  in 
Equation  (B.2). 

(x)  can  be  shown  to  be  an  unbiased  estimator  of  R^  (x)  [7],  Assuming  that  the 
process  X(t)  is  a  Gauss-Markov  process,  then  the  autocorrelation  function  is 


R,(x)=aV^'^'  (I 

Note  that  P  =  ^  •  This  autocorrelation  can  be  substituted  in  Equation  (B.l)  to  give 
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Var[V^(x)] 


(B.4) 


which  is  the  limit  that  is  used  for  the  accuracy  on  an  experimental  autocorrelation 
function  in  Chapter  3. 
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